CD Schooolhouse Version 10.0/CD Schooolhouse Version 10.0.ISO/pc/dos/math/fsultra1/ultrxcod.inc

Text File | 1992-06-18 | 8.6 KB | 375 lines

comment !
FSU - ULTRA    The greatest random number generator that ever was
        or ever will be.  Way beyond Super-Duper.
        (Just kidding, but we think its a good one.)

Authors:    Arif Zaman ([e-mail]) and
        George Marsaglia ([e-mail]).

Date:        27 May 1992

Version:    1.05

Copyright:    To obtain permission to incorporate this program into
        any commercial product, please contact the authors at
        the e-mail address given above or at

        Department of Statistics and
        Supercomputer Computations Research Institute
        Florida State University
        Tallahassee, FL 32306.

See Also:    README        for a brief description
        ULTRA.DOC    for a detailed description

-----------------------------------------------------------------------
!
.386
; This file is included in after a header file which defines
; language dependent elements. (see ultra.doc for details)
;
rinit    proc        ; INITIALIZE SEED ARRAY ========================
    Enter2arg
;
; On entry:
;    eax and ebx contain congruential and shift-register seeds
;
; On exit:
;    The swbb array is set using the McGill Super-Duper generator,
;    which is a mix of a congruential and shift-register sequence.
;    Flags and counters are reset.
;
; Registers Clobbered:  AX,BX,CX,DX,SI,DI.
;
; Algorithm:
;    Starting from lsm of x[0] to the msb of x[N-1], each bit is
;    formed using the sign bit of the xor of a congruential generator
;    with seed ConX and a shift register sequence with seed ShrX.
;
; Register usage
;   eax contains congruential seed
;   ebx contains the shift-register seed
;   cx contains count for two loops
;      ch counts the outer loop (for each byte of the array x)
;      cl counts for each bit of x[i]
;   es:di point to where x[i] is
;   edx is a temporary register
;
    mov ch,N*4
    mov di,offset swbseed    ; do loop for x[0], ... , x[4*n] (bytes)
nextbyte:
      mov cl,8
nextbit:
    mov  edx,dword ptr 69069
    mul  edx

    mov  edx,ebx
    shr  edx,15
    xor  ebx,edx
    mov  edx,ebx
    shl  edx,17
    xor  ebx,edx

    mov  edx,eax
    xor  edx,ebx
    rcl  edx,1

    rcr  byte ptr [di],1   ; shift it into the answer
    dec  cl
    jnz  nextbit
    inc  di             ; Store the answer in swbb[i]
    dec ch
    jnz nextbyte
;
; reset all counters, flags etc. and return.
;
    xor bx,bx
    mov swb32.c,bx
    mov swb16.c,bx
    mov swb8.c,bx
    mov swb1.c,bx
    mov flags,bl
    mov congx,eax
    Exit2arg
rinit    endp

SWBfill    proc near    ; SUBTRACT-WITH-BORROW ========================
;
;  On Entry:
;     DS should point to the data segment.
;     BX points to the array where the results will be stored (swb??.c)
;
;  On Exit:
;     The swbseed array contains all new values computed by the
;     subtract-with-carry generator. The CARRY byte contains the
;     state of the carry flag due to the last subtraction.
;     swb??.x contians the seed xored with a congurential.
;
;  Registers Clobbered: AX,BX,CX,DX,SI,DI
;
;  Algorithm:
;     The following subtractions are performed from right to left.
;     The carry propagates as though these were two long numbers,
;     with each x[i] being a `digit' in base 2^32.
;
;        x[12] ...  x[ 0]      x[36] ...  x[13]
;       -x[36] ... -x[24]     -x[23] ... -x[ 0]
;       ---------------------------------------
;        x[36] ...  x[24]      x[23] ...  x[ 0]
;
;
;    for a bigger table, it could be done in three steps by:
;
;     x[ 12] ...  x[   0]  :  x[N-1 ]  x[N-24]  :   x[N-25] ...  x[13]
;    -x[N-1] ... -x[N-13]  :  x[N-14] -x[N-37]  :  -x[N-38] ... -x[ 0]
;    ----------------------:--------------------:---------------------
;     x[ 36] ...  x[  24]  :  x[13  ]  x[0   ]  :   x[N-1 ] ...  x[37]
;
;     The x's also could be considered as pairs of base 16 digits,
;     so that x[i] is the pair y[2i+1]y[2i]. This allows us to use
;     only 16 bit subtractions with carry, perfectly suited for all
;     80x86 processors. The same idea could be extended for machines
;     with only eight bit, or even only 1 bit arithmetic.
;
    EnterFill
    mov ah,byte ptr flags   ; set carry flag to what it was the
    sahf                    ; last time we exited this routine

    mov cx,24               ; will do first loop 24 times
    mov ax,ds
    mov es,ax
    mov si,offset(swbseed)+13*4; set up ds:si -> x[13]
    mov di,offset(swbseed)     ; set up es:di -> x[0]

loop1:  ; On a 80386, the instructions `lodsd', `stosd' and `loop'
    ; all change registers automatically as noted in parentheses

    lodsd                   ;    ax = x[i+13]     ( si = si+4 )
    sbb eax,[di]            ;    ax = ax-x[i]-carry
    stosd                   ;    y[i] = ax        ( di = di+4 )
    loop loop1              ; loop for i=0..47    ( cx = cx-1 )

    mov cx,13               ; will do next loop 13 times
    mov si,offset(swbseed)  ; set up ds:si -> x[0]
                                ; es:di is already set up
loop2:
    lodsd                   ;    ax = x[i-24]     ( si = si+4 )
    sbb eax,[di]            ;    ax = ax-x[i]-carry
    stosd                   ;    x[i] = ax        ( di = di+4 )
    loop loop2              ; loop for i=48..73   ( cx = cx-1 )

    lahf
    mov byte ptr flags,ah   ; save carry flag for next time
;
; XOR the elements of swbb with a congruential generator and put the
; result in swbnn.x, reset the counter and the pointer.
;

    mov  di,bx
    mov  [bx-4],bx
    mov  si,offset(swbseed)
    mov  cx,N
IFNDEF fast
    mov  eax,congx
    mov  ebx,69069
loopc:  mul  ebx
    mov  edx,eax
    lodsd
    xor  eax,edx
    stosd
    mov  eax,edx
    loop loopc
    mov  congx,eax
ELSE
    rep  movsd
ENDIF
    ExitFill
SWBfill    endp

; Random Number procedures ============================================
;
CheckFill MACRO bits,bytes,count
    local ok
    dec    bits.c
    jns    ok
        mov   bx,offset(bits.x)
        call  SWBfill
        mov   bits.c,count-1
ok:    mov    bx,bits.p
    add    bits.p,bytes
ENDM

i32bit proc
    EnterProcedure
    CheckFill swb32,4,N
    mov    ax,[bx]
    mov    dx,[bx+2]
    DwordFn
    ExitProcedure
i32bit endp

i31bit proc
    EnterProcedure
    CheckFill  swb32,4,N
    mov    ax,[bx]
    mov    dx,[bx+2]
    and    dh,7Fh
    DwordFn
    ExitProcedure
i31bit endp

i16bit proc
    EnterProcedure
    CheckFill  swb16,2,2*N
    mov    ax,[bx]
    WordFn
    ExitProcedure
i16bit endp

i15bit proc
    EnterProcedure
    CheckFill  swb16,2,2*N
    mov    ax,[bx]
    and  ah,7Fh
    WordFn
    ExitProcedure
i15bit endp

i8bit  proc
    EnterProcedure
    CheckFill  swb8,1,4*N
    mov    al,[bx]
    ByteFn
    ExitProcedure
i8bit  endp

i7bit proc
    EnterProcedure
    CheckFill  swb8,1,4*N
    mov    al,[bx]
    and al,7Fh
    ByteFn
    ExitProcedure
i7bit  endp

i1bit  proc
    EnterProcedure di
    dec    swb1.c
    jns    ok1
        CheckFill  swb32,4,N    ; do an i32bit
        mov      dx,[bx+2]
        mov   bx,[bx]        ; with the answer in dx:bx
        mov   swb1.c,31
        mov   di,offset(swb1.x)
        mov   swb1.p,di
        mov   ax,ds
        mov   es,ax

        mov   cx,16
stosb1:        mov      al,1
        and   al,bl
        stosb
        shr      bx,1
        loop  stosb1

        mov      cl,16
stosb2:        mov   al,1
        and   al,dl
        stosb
        shr   dx,1
        loop  stosb2
ok1:    mov    bx,swb1.p
    inc    swb1.p
    mov    al,[bx]
    ByteFn
    ExitProcedure di
i1bit  endp

uni proc
    EnterProcedure
    fild    neg31
    CheckFill    swb32,4,N
    and    byte ptr [bx+3],7Fh
    jnz    oku
        mov   eax,[bx]
        mov   tmpdhi,eax
        CheckFill  swb32,4,N
        mov   eax,[bx]
        mov   tmpdlo,eax
        fild  tmpq
        jmp   uxit
oku:    fild    dword ptr [bx]
uxit:    fscale
    fstp    st(1)
    RealFn
    ExitProcedure
uni endp

vni proc
    EnterProcedure
    fild    neg31
    CheckFill    swb32,4,N
    test    byte ptr [bx+3],0FFh
    jnz    okv
        mov   eax,[bx]
        mov   tmpdhi,eax
        CheckFill  swb32,4,N
        mov   eax,[bx]
        mov   tmpdlo,eax
        fild  tmpq
        jmp   vxit
okv:    fild    dword ptr [bx]
vxit:    fscale
    fstp    st(1)
    RealFn
    ExitProcedure
vni endp

duni proc
    EnterProcedure
    fild    neg63
    CheckFill    swb32,4,N
    dec    swb32.c
    jns    okdu
        mov   eax,swb32.x[4*N-4]
        mov   tmpdlo,eax
        mov   bx,offset(swb32.x)
        call  SWBfill
        mov   swb32.c,N-1
        mov      eax,swb32.x
        mov   tmpdhi,eax
        and   byte ptr [tmpq+7],7Fh
        fild  tmpq
        jmp   duxit
okdu:    and    byte ptr [bx+7],7Fh
    fild    qword ptr [bx]
duxit:    add    swb32.p,4
    fscale
    fstp    st(1)
    DoubleFn
    ExitProcedure
duni endp

dvni proc
    EnterProcedure
    fild    neg63
    CheckFill  swb32,4,N
    dec    swb32.c
    jns    okdv
        mov   eax,swb32.x[4*N-4]
        mov   tmpdlo,eax
        mov   bx,offset(swb32.x)
        call  SWBfill
        mov   swb32.c,N-1
        mov   eax,swb32.x
        mov   tmpdhi,eax
        fild  tmpq
        jmp   dvxit
okdv:    fild    qword ptr [bx]
dvxit:    add    swb32.p,4
    fscale
    fstp    st(1)
    DoubleFn
    ExitProcedure
dvni endp